home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / bessel_I1.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  6.5 KB  |  259 lines
  1. /* specfunc/bessel_I1.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_bessel.h>
  26.  
  27. #include "error.h"
  28.  
  29. #include "chebyshev.h"
  30. #include "cheb_eval.c"
  31.  
  32. #define ROOT_EIGHT (2.0*M_SQRT2)
  33.  
  34.  
  35. /*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/
  36.  
  37. /* based on SLATEC besi1(), besi1e() */
  38.  
  39. /* chebyshev expansions
  40.  
  41.  series for bi1        on the interval  0.        to  9.00000d+00
  42.                     with weighted error   2.40e-17
  43.                      log weighted error  16.62
  44.                    significant figures required  16.23
  45.                     decimal places required  17.14
  46.  
  47.  series for ai1        on the interval  1.25000d-01 to  3.33333d-01
  48.                     with weighted error   6.98e-17
  49.                      log weighted error  16.16
  50.                    significant figures required  14.53
  51.                     decimal places required  16.82
  52.  
  53.  series for ai12       on the interval  0.        to  1.25000d-01
  54.                        with weighted error   3.55e-17
  55.                     log weighted error  16.45
  56.                   significant figures required  14.69
  57.                    decimal places required  17.12
  58. */
  59.  
  60. static double bi1_data[11] = {
  61.   -0.001971713261099859,
  62.    0.407348876675464810,
  63.    0.034838994299959456,
  64.    0.001545394556300123,
  65.    0.000041888521098377,
  66.    0.000000764902676483,
  67.    0.000000010042493924,
  68.    0.000000000099322077,
  69.    0.000000000000766380,
  70.    0.000000000000004741,
  71.    0.000000000000000024
  72. };
  73. static cheb_series bi1_cs = {
  74.   bi1_data,
  75.   10,
  76.   -1, 1,
  77.   10
  78. };
  79.  
  80. static double ai1_data[21] = {
  81.   -0.02846744181881479,
  82.   -0.01922953231443221,
  83.   -0.00061151858579437,
  84.   -0.00002069971253350,
  85.    0.00000858561914581,
  86.    0.00000104949824671,
  87.   -0.00000029183389184,
  88.   -0.00000001559378146,
  89.    0.00000001318012367,
  90.   -0.00000000144842341,
  91.   -0.00000000029085122,
  92.    0.00000000012663889,
  93.   -0.00000000001664947,
  94.   -0.00000000000166665,
  95.    0.00000000000124260,
  96.   -0.00000000000027315,
  97.    0.00000000000002023,
  98.    0.00000000000000730,
  99.   -0.00000000000000333,
  100.    0.00000000000000071,
  101.   -0.00000000000000006
  102. };
  103. static cheb_series ai1_cs = {
  104.   ai1_data,
  105.   20,
  106.   -1, 1,
  107.   11
  108. };
  109.  
  110. static double ai12_data[22] = {
  111.    0.02857623501828014,
  112.   -0.00976109749136147,
  113.   -0.00011058893876263,
  114.   -0.00000388256480887,
  115.   -0.00000025122362377,
  116.   -0.00000002631468847,
  117.   -0.00000000383538039,
  118.   -0.00000000055897433,
  119.   -0.00000000001897495,
  120.    0.00000000003252602,
  121.    0.00000000001412580,
  122.    0.00000000000203564,
  123.   -0.00000000000071985,
  124.   -0.00000000000040836,
  125.   -0.00000000000002101,
  126.    0.00000000000004273,
  127.    0.00000000000001041,
  128.   -0.00000000000000382,
  129.   -0.00000000000000186,
  130.    0.00000000000000033,
  131.    0.00000000000000028,
  132.   -0.00000000000000003
  133. };
  134. static cheb_series ai12_cs = {
  135.   ai12_data,
  136.   21,
  137.   -1, 1,
  138.   9
  139. };
  140.  
  141.  
  142. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  143.  
  144. int gsl_sf_bessel_I1_scaled_e(const double x, gsl_sf_result * result)
  145. {
  146.   const double xmin    = 2.0 * GSL_DBL_MIN;
  147.   const double x_small = ROOT_EIGHT * GSL_SQRT_DBL_EPSILON;
  148.   const double y = fabs(x);
  149.  
  150.   /* CHECK_POINTER(result) */
  151.  
  152.   if(y == 0.0) {
  153.     result->val = 0.0;
  154.     result->err = 0.0;
  155.     return GSL_SUCCESS;
  156.   }
  157.   else if(y < xmin) {
  158.     UNDERFLOW_ERROR(result);
  159.   }
  160.   else if(y < x_small) {
  161.     result->val = 0.5*x;
  162.     result->err = 0.0;
  163.     return GSL_SUCCESS;
  164.   }
  165.   else if(y <= 3.0) {
  166.     const double ey = exp(-y);
  167.     gsl_sf_result c;
  168.     cheb_eval_e(&bi1_cs, y*y/4.5-1.0, &c);
  169.     result->val  = x * ey * (0.875 + c.val);
  170.     result->err  = ey * c.err + y * GSL_DBL_EPSILON * fabs(result->val);
  171.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  172.     return GSL_SUCCESS;
  173.   }
  174.   else if(y <= 8.0) {
  175.     const double sy = sqrt(y);
  176.     gsl_sf_result c;
  177.     double b;
  178.     double s;
  179.     cheb_eval_e(&ai1_cs, (48.0/y-11.0)/5.0, &c);
  180.     b = (0.375 + c.val) / sy;
  181.     s = (x > 0.0 ? 1.0 : -1.0);
  182.     result->val  = s * b;
  183.     result->err  = c.err / sy;
  184.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  185.     return GSL_SUCCESS;
  186.   }
  187.   else {
  188.     const double sy = sqrt(y);
  189.     gsl_sf_result c;
  190.     double b;
  191.     double s;
  192.     cheb_eval_e(&ai12_cs, 16.0/y-1.0, &c);
  193.     b = (0.375 + c.val) / sy;
  194.     s = (x > 0.0 ? 1.0 : -1.0);
  195.     result->val  = s * b;
  196.     result->err  = c.err / sy;
  197.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  198.     return GSL_SUCCESS;
  199.   }
  200. }
  201.  
  202.  
  203. int gsl_sf_bessel_I1_e(const double x, gsl_sf_result * result)
  204. {
  205.   const double xmin    = 2.0 * GSL_DBL_MIN;
  206.   const double x_small = ROOT_EIGHT * GSL_SQRT_DBL_EPSILON;
  207.   const double y = fabs(x);
  208.  
  209.   /* CHECK_POINTER(result) */
  210.  
  211.   if(y == 0.0) {
  212.     result->val = 0.0;
  213.     result->err = 0.0;
  214.     return GSL_SUCCESS;
  215.   }
  216.   else if(y < xmin) {
  217.     UNDERFLOW_ERROR(result);
  218.   }
  219.   else if(y < x_small) {
  220.     result->val = 0.5*x;
  221.     result->err = 0.0;
  222.     return GSL_SUCCESS;
  223.   }
  224.   else if(y <= 3.0) {
  225.     gsl_sf_result c;
  226.     cheb_eval_e(&bi1_cs, y*y/4.5-1.0, &c);
  227.     result->val  = x * (0.875 + c.val);
  228.     result->err  = y * c.err;
  229.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  230.     return GSL_SUCCESS;
  231.   }
  232.   else if(y < GSL_LOG_DBL_MAX) {
  233.     const double ey = exp(y);
  234.     gsl_sf_result I1_scaled;
  235.     gsl_sf_bessel_I1_scaled_e(x, &I1_scaled);
  236.     result->val  = ey * I1_scaled.val;
  237.     result->err  = ey * I1_scaled.err + y * GSL_DBL_EPSILON * fabs(result->val);
  238.     result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
  239.     return GSL_SUCCESS;
  240.   }
  241.   else {
  242.     OVERFLOW_ERROR(result);
  243.   }
  244. }
  245.  
  246. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  247.  
  248. #include "eval.h"
  249.  
  250. double gsl_sf_bessel_I1_scaled(const double x)
  251. {
  252.   EVAL_RESULT(gsl_sf_bessel_I1_scaled_e(x, &result));
  253. }
  254.  
  255. double gsl_sf_bessel_I1(const double x)
  256. {
  257.   EVAL_RESULT(gsl_sf_bessel_I1_e(x, &result));
  258. }
  259.